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ABSTRACT 

The term "vortex breakdown" refers to the abrupt and drastic changes of 
structure that can sometimes occur in swirling flows. It has been conjectured that 
the "bubble" type of breakdown can be viewed as an axisymmetric wave travel- 
ling upstream in a primarily columnar vortex flow. In this scenario the wave’s 
upstream progress is impeded only when it reaches a critical amplitude and it 
loses stability to some non- axisymmetric disturbance. We will investigate the sta- 
bility of some axisymmetric wavy flows, which model vortex breakdown, to three 
dimensional disturbances viewing the amplitude of the wave as a bifurcation 
parameter. We will also look at the stability of a set of related, columnar vortex 
flows which eire constructed by taking the two dimensional flow at a single axial 
location and extending it throughout the domain without variation. The method 
of our investigation will be to expand the perturbation velocity in a series of diver- 
gence free vectors which ensures that the continuity equation for the incompressi- 
ble fluid is satisfied exactly by the computed velocity field. Projections of the sta- 
bility equation onto the space of inviscid vector fields eliminates the pressure term 
from the equation and reduces the differential eigen problem to a generalized 
matrix eigen problem. Results are presented both for the one dimensional, colum- 
nar vortex flows and also for the wavy "bubble" flows. 
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1. Introduction: Vortex Breakdown. 


The term "vortex breakdown" refers to the abrupt and drastic changes of structure that can 
sometimes occur in vortex flows. Observations by Peckham & Atkinson [1957] of breakdowns 
occurring in the leading edge vortex formed above a swept back lifting surface and a number of 
studies demonstrating the serious aerodynamic consequences of such events (the slopes of the lift, 
drag and moment curves are all altered by breakdown) stimulated early interest in the subject. 
Since that time the literature on vortex breakdown has burgeoned. The interested reader is 
referred to review articles by Hall [1972] and Leibovich [1978, 1984] for summaries both of the 
experimental observations that have been made and also of the theories that have been proposed 
to explain them. 

Experimental observations are most easily made on vortex flows confined to tubes and the 
bulk of the available data is for such cases. In one apparatus, used by a number of researchers, 
water is passed radially inward through a set of guidevanes imparting swirl to the fluid which 
then enters axially into a test section (a frustrum of a cone of very small cone angle), by means of 
an annular channel formed between a bellmouth opening on the section and a centerbody whose 
tip is aligned with the cone axis. The boundary layer shed from the tip of the centerbody forms a 
well defined vortex core along the axis of the test section and dye injected through the tip allows 
for flow visualization. 

With this type af apparatus two parameters are within the easy control of the experimen- 
talist, namely the amount of swirl imparted to the inlet flow and the volume flow rate through 
the tube (effectively the Reynolds number of the flow). As the Reynolds number is increased for 
a sufficently large, fixed value of swirl the breakdown assumes one of two characteristic forms. 
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Both of these are characterized by a rapid deceleration of the axial velocity component, occurring 
in the axial distance on the order of one vortex core diameter, followed by the formation of a 
stagnation point and (in some frame of reference) a region of reversed flow along the axis. The 
two forms are easily distinguished in flow visualization studies as in one form, the spiral or S type 
breakdown, the tracer dye assumes a spiral shape rotating in the same sense as the inlet fluid, 
while in the other form, the bubble or B type breakdown, the dye assumes a form that looks 
much like a body of revolution placed in the fluid. Our interest will be in this latter form of 
breakdown which is sketched in Figure 1. Here we show "ideal" or averaged stream surfaces on 
which the fluid particles travel in helical paths. (Leibovich [1978]). 

Faler &. Leibovich [1977], Garg & Leibovich [1979] and the author [unpublished studies] 
have used the non-invasive techniques of laser doppler anemometry to measure the velocity fields 
both upstream and downstream of breakdown events. Outside a thin boundary layer along the 
tube wall the experimental data is well fitted by the analytic profiles, 

F(r) - ifl(l - <T"') (1.1) 

W(r) = W 1 + W 2 e~ ar * (1.2) 

W and V being respectively the axial and azimuthal velocity components while W 2 , Q and a 
are all constants (representative values are given by Garg & Leibovich [1979]). The profiles apply 
to the downstream flow only in the mean, as the flow there fluctuates with time. 

Upstream of the recirculation zone the flow is axisymmetric and steady. After breakdown 
the vortex core expands to two or three times its upstream size and the constant W 2 in the mean 
axial velocity profile which had been positive upstream (jetlike flow) becomes negative (wakelike 
flow). Downstream, within a few vortex core diameters of the breakdown, a turbulent wake is 
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in variably established - This transition to turbulence "switched on" by the coherent breakdown 
structure provides a further incentive for its study. 

Possibly motivated by the fact that the flows upstream of breakdown can be made to have a 
high degree of axial symmetry, most of the research to date assumes that axially symmetric 
processes are the important ones in vortex breakdown. As only axisymmetric disturbances can 
cause a change in the axial velocity component as measured on the axis and as a deceleration of 
this component is so pronounced in breakdown, it is clear that such disturbances play an impor- 
tant role. Nevertheless, all transitions occurring in vortex flows as documented by Faler [1976] 
are nonaxisymmetric and the flow within the bubble itself is unsteady with regular low frequency 
oscillations. Furthermore, the stagnation point that defines the start of the recirculation zone is 
not entirely fixed but wanders over a short range of the axis in a seemingly random fashion. 

Leibovich [1984] proposed the following plausible scenario for the bubble type breakdown. 
A finite axisymmetric disturbance, triggered off downstream, moves upstream in a columnar flow 
that is nearly critical in the sense of Benjamin [1962]. (A supercritical flow, in this classification, 
allows for the upstream propagation of infinitesimal axisymmetric waves while a subcritical flow 
does not). Flows of the form (1.1,2) can indeed support axisymmetric dispersive waves and these 
can propagate upstream in some situations (Leibovich [1970], Randall & Leibovich [1973]). Mov- 
ing in this direction, the cross sectional area of the tube decreases causing the wave to amplify 
and speed up. The conjecture is that, upon reaching some critical amplitude, these waves lose 
stability to a non-axisymmetric disturbance. The growth of the asymmetric mode at the expense 
of the axisymmetric wave, drains energy from it and this causes the wave to equilibrate at some 
axial location in the diverging tube, much as is seen in experiments. 
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Our aim is to study the stability of some inviscid, wavy axisymmetric flows to three dimen- 
sional disturbances with the amplitude of the waves viewed as a bifurcation parameter. We start 
with a columnar flow in cylinderical coordinates of the form (0, V 0 (r), lV 0 (r)), ( e -g- (1-1,2)). (For 
arbitrary C l functions, V 0 and W 0 , all such flows satisfy Euler’s equations). We then seek 
axisymmetric wavy perturbations to this flow which satisfy the equations of motion, at least 
approximately, for small amplitude. In terms of a stream function, tp and a circulation, k 
(Lamb [1932]) Leibovich [1972] found solutions to Eulers equations of the form, 



= M r ) + ^4>{r)A{x) + 0(e 2 ), 

(1.3) 

n(r,x) 

= K 0 (r) + e 7 (r)yi(x), + 0(e 2 ), 

(1.4) 

where x is a moving coordinate, 




x = z — dt 

(1.5) 

and d is a constant that must found in the calculation. The velocity components 

are given by 


1 9 , 

u = — tp, 

r az 

(1.6) 


1 

V = K, 

r 

(1.7) 


i a , 

w = — — 1 />. 

r dr 

(1.8) 


The columnar base flow is represented by V’o( r ) an ^ K o[ r )‘ 

The amplitude function, A(z,t) is governed by a Korteweg de Vries equation which has both 
infinite and finite period solutions. The multiple scales analysis that was used to obtain these 
solutions is strictly valid only for long period waves which Eire also the most interesting solutions 
from a physical point of view. When doing the stability analysis we will confine our attention, 
for numerical reasons, to solutions of the finite period, 2 H and these are given exactly in terms of 
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cnoidal functions (Whitham [1974]). The structure functions, 4>(r), 7 (r) and the wave speed d 
are determined (numerically) from a second order, ordinary differential eigenvalue problem. 

For certain base columnar profiles d is negative and consequently the axisymetric wave pro- 
* pagates upstream. Figure 2 is a plot of the streamlines (1.3) in a frame moving with the wave for 

such a case. The base columnar profile used here and throughout this paper is a purely swirling 
flow; W 0 (r) = 0 and a - 14 in the notation of (1.1). The structure function <f> has been normal- 
ized so that Max <f> = 1 and for this flow a recirculation zone (bubble) first appears in the stream- 
line plot when the amplitude parameter, e reaches a value of 0.0155. For the value of e used here 
the plot is clearly reminiscent of the bubble type breakdown. 

Our aim is to study the stability of the flows (1.3,4) to three dimensional disturbances view- 
ing e as a bifurcation parameter. The analysis will be carried out in a frame moving with the 
wave, i.e. using the coordinates (r, x, 6). As the base flow is dependent on both the radial and 
axial variables, r and 2 , the stability equations separate only in the azimuthal variable, 9. It will 
be in our interest also to study the stability of a related columnar flow that is constructed by tak- 
ing the two dimensional flow (1.3,4) at a single axial station, x = 0, and extending it throughout 
the cylindrical domain without variation. For obvious reasons we will refer to this flow as the 
"mid-bubble" columnar flow and it is given explicitly as follows, 


A 

nw - 

w 

+ 77 ( 0 , 

(1.9) 

jL 

w.M - 

W 0 (r) 

+ -V(r). 

(1.10) 


Plots of these profiles for various values of e are given in Figures 3 and 4. Provided the wavy 
flow (1.3,4) varies only slowly along the axis (as it will do if the period, 2 H of A(x) is very large), 
we can look on the midbubble profiles as models for the full two dimensional flow. We conjecture 
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that the stability of these columnar flows (the equations for which separate in both r and 0 ) 
should also be indicative of the stability properties of the full two dimensional flow. 

In the rest of this paper we describe the numerical scheme used to solve the stability equa- 
tions, we discuss its implementation and verification on the computer and finally we give results 
obtained for the stability of the midbubble columnar and the axisymmetric wavy flows presented 
above. 


2. Numerical Methods. 

For incompressible fluids the physical law of mass conservation reduces to the constraint that the 
velocity vector of the fluid be divergence free. The pressure is not then a thermodynamic variable 
determined by an equation of state but rather can be thought of as a Lagrange multiplier adjust- 
ing itself instantaneously to ensure that this kinematical constraint on the velocity vector is met. 
There is no evolution equation for the pressure nor does it satisfy any predetermined boundary or 
initial conditions. 

Numericists, seeking to solve the governing equations approximately, have found that their 
greatest difficulty lies in the treatment of the pressure variable. While many ingenious methods 
have been devised to overcome the difficulties, the treatment advocated in this work is in a 
mathematical sense the most natural and offers many computational advantages. Here, the pres- 
sure term is eliminated from the equations entirely and the divergence free condition is satisfied 
exactly by the numerically obtained approximation to the velocity vector. Moreover, as the com- 
ponents of the velocity are expanded in terms of series of polynomials that arise as the solution to 
a singular Sturm Liouville problem, whose excellent approximation properties are well 
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documented (e.g. Gottleib & Orszag [1977], Quarteroni [1983]) convergence of our approximation 
will be bound only by the smoothness of the solution and by the number of terms used in the 
component expansions. For infinitely differentiable velocity fields we should expect to achieve 
"exponential convergence" (Canuto et al.) 

The essence of the method (originally due to Leonard &; Wray [1982]) involves expanding 
the velocity in a series of divergence free vector fields each of which satisfy the same boundary 
conditions as the velocity. The infinite sums are truncated and substituted into the governing 
equations, which are the Navier-Stokes or Euler equations linearized about the appropriate base 
flow. Inner products are taken with vectors fields which satisfy inviscid boundary conditions. 
This eliminates the pressure term from the equations and reduces the differential eigenvalue prob- 
lem to a matrix eigenvalue problem. The eigenvalues determine the stability of the base flow and 
the eigenfunctions are the set of coefficents in the expansions of the corresponding perturbation 
velocity fields. 

To examine how this method works we recall that it is well known (Ladyshenskaya [1966]) 
that L 2 {D), the space of square integrable vector functions defined on a bounded domain 
D{Dc R n ’ n = 2,3) can be decomposed into those that are divergence free and whose normal 
components vanish on the boundary and those that can be expressed as the gradient of a 
differentiable function defined on D. For this paper we will consider vector fields, defined on the 
section of a cylinder T, which are periodic in both the axial and azimuthal variables, having as 
their axial period the tube length, 2H. 

We will decompose L 2 {T) as follows. 

L\T) = J{T) + J l {T), 


where, 


( 2 . 1 ) 
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J{T) = 


(a) V‘u = 0 in T, 
u e L 2 (T) (b) u-n = 0 on 3T, 

( c ) “i 5 , = «l 5 , • 


( 2 . 2 ) 


S 1 and S 2 represent the ends of the cylinder. Given in this form J(T) is clearly the space of 
(a) incompressible, (b) inviscid, (c) periodic velocity fields. 

The set of "viscous" velocity fields on T is a subset of J(T) denoted J°(T). 

J°[T) = jueJ(r) | u = 0on3r|. (2.3) 

An alternative representation of J(T ) (Richtmyer [1978]) is given by, 


J(T) = 


ueL 2 {T) 


(a) <u, VP> = 0 f or a U P £ C°°(T) 


( 6 ) «1 5, = «I 5, 

where T is the closure of T and <*,•> represents the usual inner product in L 2 (T), 


(2.4) 


<“»«> = It u*£ rdrdOdx . (2-5) 

The space J(T ) endowed with this inner product is a Hilbert space and a closed subspace of 
L 2 (T). The projection of L 2 (T ) onto J(T) will be denoted by II. It is clear that vectors of the 
form VP are perpendicular to all u in J(T) and in fact (Ladyshenskaya [1966]), 


J^(T) = | u e L 2 (T) | u = Vp for some p in C^T) | . 


( 2 . 6 ) 


II then has the following properties, 


n : L\T) -> J{T), 


(2.7) 
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( 2 . 8 ) 

(2.9) 


II u = u for all « c J(T) , 

n yp = 0 for all P e C\T). 

To determine the linear stability of a flow U to say, viscous disturbances which are periodic 
in x and 6 we consider whether infinitesimal perturbations to U grow in time. Therefore we 
linearize the Navier-Stokes equations about U and seek solutions in J°(T) of the form, 

u(r,x,0)e~ ,ffl . (2.10) 

The character of <r then determines the linear temporal stability of U. If a = a + :/? where a, /? 
are real then, 

/? > 0 => U is unstable, 

P = 0 => tP is neutrally stable,. (2. 11) 

f3 < 0 => U is stable 

The equations that must be solved have the form, 

iVu = Eu + Re -1 5u. (2.12) 

E and S are operators defined on J[T) as follows, 

Su = -n(yx^) (2.13) 

and 

Eu = II(wxCl + Hxu), (2.14) 

where fi and w are respectively the base and perturbation vorticities, (fl = Vx t/, uj_ = yxti). 

Some suitable nondimensionalization has introduced the Reynolds number, 

UnRn 

R e =_L_2. j (2.15) 

v 

R q and U 0 being chaxacteristic length and velocity scales associated with the base flow and u is 
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the kinematic viscosity of the fluid. We can take R 0 to be the radius of the tube and U 0 to be the 
maximum value of the columnar base flow azimuthal velocity, Vo(r). 

The VP term in the Navier-Stokes equation has been eliminated by projection onto J(T). 
Projection of the stability equation onto a finite dimensional subspace, of J(T) is achieved 

in practice by taking the inner product of the equation with basis vectors for J N (T). This pro- 
cess eliminates the operator II from the equation, for given any vector / in L 2 (T) and any vector 
A in J(T) we have that, 


<11 /, A> = </, A>, (2-16) 

as projections are self adjoint and as the projection of any vector in J(T) is itself. 

It is worth emphasising that even when we are solving the viscous stability equations, we 
still project the governing equations onto the space of inviscid vector fields. The reason for this is 
that having found a velocity u such that the vector / defined as, 

/ = icru — wxP — nxu + j/V x |£. (2-17) 

is orthogonal to all A in J(T) then / f J^(T) and so there exists a scalar function p (a pressure) 
with / = vp. If? however, / were in which contains J^(T) then the existence of a pres- 

sure is not guaranteed and consequently u may not correspond to a physical solution. 

Leonard and Wray [1982] demonstrated a divergence free vector function expansion for 
viscous velocity fields, defined on a cylindrical domain that are Fourier decomposable in both the 
axial and azimuthal variables. We will construct a somewhat different set of basis vectors here. 
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The velocity field, u, satisfies the continuity equation and is Fourier decomposable in x and 
6 , which means in effect that only two of its three components, u, v, w are independent. This 
motivates the introduction of two vector families, in an expansion of the form, 


« = + fl 2n-UmX»( r )| c + 

nkm \ ) 


(2.18) 


The components of the vectors are found as follows. Expand two of the velocity components, 
say the first and the third, independently as, 


nkm 


nkm 


(2.19) 


( 2 . 20 ) 


where f„{r) are complete sets of polynomials chosen to satisfy the boundary conditions that are 
imposed on tt, w. The r and x components of xjf have now been picked and it remains for us to 
chose the 6 components in a manner that ensures the vectors x % e + are divergence free. 
Consider for example, xJT(r). 


=> 


v{x;(r)e*^ + ^)j = 0, 


( r /»( r ))' + im X n ,e = 0, 


( 2 . 21 ) 


( 2 . 22 ) 


where the prime denotes a derivative with respect to r. This equation gives us the 6 component of 
XZ( r )' Rescaling, it is found that an expansion of the form (2.14) is possible for non-zero azimu- 
thal wavenumbers where, 




(r) = | i'm/„ (»•)> ~ {rf n ( r )) '» °) > 


(2.23) 
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ti(r) = (o,-r*/ B + (r),m/ B + (r)) (2.24) 

and such an expansion will guarantee that u is divergence free. This expansion is clearly incom- 
plete for azimuthal wave number zero, (tn = 0), i.e. for axisymmetric flows. For that case the 
following expansion vectors can be used. 


X»( r ) = 


**/„(»•)> 0, - — (r/ B (r))' 


(2.25) 


- (o, /. + (r),o) (2.26) 

The polynomials / B (r) must be chosen so that the vector u given by (2.14) satisfies 
appropriate (viscous or inviscid) boundary conditions on the walls of the domain, T and is single 
valued at the origin, r = 0. We will denote the polynomials used in the inviscid case by a B (r) 
reserving / B (r) for viscous expansions. We have then upon truncating (2.14) an approximation 
to u of the form, 

N K M 

htfKM = S 2 S a nkm ^n*m( r 5 2 )^) (2.27) 

n=l k-—Km=—M 

where, 

£*»(»•,*, *) * X^(r-,k,m) e < ^ + me l (2.28) 

The projection vectors will have the same form as the expansion vectors, i.e. we will project 
with vectors, Aj pf (r,z,0) , where 


for 


A Jf J r ,x,0) = |f(r + 


I = l,...,jV; p = q = 


(2.29) 


with the vectors being given by equations (2.23 - 26) using the inviscid polynomials, o*(r) for 
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the components. 

It is possible to choose the polynomials a*(r) and /*(r) in many different ways. Leonard & 
Wray [1982] in their consideration of certain turbulence simulations employed an unusual set of 
Jacobi polynomials to reduce the bandwidth of the final matrix system. These polynomials were 
also used by Spalart [1983] in his simulation of boundary-layer transition. Moser & Moin [1984] 
in their work on the infinite Taylor Couette system, used Tchebychev polynomials and incor- 
porated the weight function, against which these polynomials are orthogonal, into the projection 
vectors. Here, we will construct the basis vectors from Legendre polynomials. All of the above 
sets are solutions to singular Sturm Liouville problems and consequently we can expect expan- 
sions in terms of any of these polynomials to exhibit excellent convergence properties. 

The single valuedness criterion, which must be applied along the centre line of the tube for 
the vector u, causes the polynomials /*(r) and af(r) to depend on m, the azimuthal 
wavenumber (Joseph [1970]). One appropriate choice for o ; ± (r) is, 

<h + (r) = rP,( 2r - 1) for all 

a f ( r ) = (! - 0P/(2r - 1) if M = 1, (2.30) 

af (r) = r(l - r)P,{ 2r - 1) if I m l * 

where the radial variable has been scaled by the tube radius and Pj(r) is the Legendre polynomial 
of order l (Abramowitz & Stegun [1970]). The corresponding choice for /*(r) is, 

fnir) = r(l - r)P n (2r - 1) for all 


/»(0 = (1 - r) 2 P n (2r - 1) if | m| = 1, 
/»(r) = r(l - r) 2 P n (2r - 1) ^M^l. 


(2.31) 
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As all of our stability problems separate in the azimuthal direction, this dependence on m 
presents no difficulty. We solve separate problems for each azimuthal wavenumber; so having 
chosen an m the expansion and projection sets are fixed throughout the calculation . Indeed, in 
theory there is no difficulty even if the problem at hand is truly three dimensional; however some 
care is required in implementing the method to ensure that the correct radial polynomial set is 
being used for each azimuthal component of the velocity. 


3. Implementation and Verification of the Method. 


Equation (2.12) is solved approximately by using the expansion f° r « and taking inner pro- 
ducts of the equation with the projection vectors, Aj pJ to get a generalized matrix eigen problem 
for the eigenvalues a and the eigenvectors a (the coefficents in the expansion u^km)- This matrix 
problem can be written as, 


crAa 


D i t 


B + 


Re 


(3.1) 


The matrix A is purely real and arises from the fact that the expansion and projection vectors 
are not orthonormal. 


Alnpkim ^^JpV Xji^^pk^jmi ~ (^- 2 ) 

The Kronecker delta symbol, S t arises because the Fourier bases employed in the axial and 
azimuthal directions are orthogonal. The matrix B arising from the convection terms is also 
purely real. 

^Inpkqm ~ < 5 (V* fiitm) * ^ d" (*^'*0 

Finally, the matrix C arising from the viscous terms is purely imaginary. 
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^Inpkqm (3-4) 

Using the orthogonality of the Fourier bases it can be written as, 

^Inpkqm = ^In^pk^qm. ( 3 - 5 ) 

The form of the matrix B depends on the base flow U. F or columnar flows which are independent 
of x and 8 it is possible to find a matrix B such that, 

^Inpkqm = Bln^pk^im. ( 3 - 6 ) 

The stability of these flows can be determined by solving the O(N) generalized matrix eigen 
problem, 


crAa = 


B + ± c 


(3.7) 


On the other hand, for the axisymmetric wavy base flow (1.3,4) we have 


^Inpkqm ^Inpk^pk 


and the stability equation is the C>| (2 K + l)xivj matrix equation, 


cAa = 


B + Re ° 


(3.8) 


(3.9) 


where A and C are the (2 K + 1) x ivj block diagonal matrices A n^pk and C^Spk respectively. 

The matrices depend parameterically on the wavenumbers of the projection and expansion vec- 
tors so we separate them into submatrices that can be evaluated independently of these and the 
other parameters (in particular e) occurring in the base flow. The submatrices are evaluated once 
and then stored in the computer. The required integrations can be done at very little cost by util- 
izing the orthogonality properties of the expansion and projection polynomials. The full matrices 
are then be reassembled without the need for doing any further integrations. 
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One can always band the A and C matrices by appropriate choice of the polynomials /*(r) 
and a/ t (r). However the matrix B will generally be full, though for certain rather simple base 
flows such as the Hagen Pouiseille flow considered be Leonard & Wray [1982] it is also possible to 
band B. The matrix A was inverted to produce a regular eigenvalue problem in place of (3.1) 
and the QR algorithm was used to extract the eigenvalues. We also note that the matrix problem 
we get when considering the inviscid stability of base flows is a purely real one and consequently 
the eigenvalues occur, as they should do, in conjugate pairs. 

A computer code has been written which implements the method we have been describing to 
solve the stability problems, both viscous and inviscid, for all columnar flows and for axisym- 
metric flows of the form (1.3,4). Both the direct and adjoint versions of the stability problems 
can be solved. The adjoint viscous stability problem is to find a u in J°[T) such that 

iXu = E*u + (3.10) 

The operator E* is the adjoint operator to E and is given by, 

e \ = -n( nxu + vx(«xt/)j (3.ii) 

The direct and adjoint spectra obtained by solving (2.12) and (3.10) should, of course, be conju- 
gate to each other and how well a numerical scheme reproduces this theoretical result is a test of 
its accuracy. 

We verified the code by calculating the stability of rotating Poiseuille flow, 

U = ( 0, V, ^i(l - r 2) ) (3.12) 

Cotton et al. [1980] found that this flow with V x = 0.2147 and W x = 1.0 was neutrally stable to 
disturbances having azimuthal wavenumber, m = 1 and axial wavenumber, k — — 1 for a 
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Reynolds number of 156. The following table lists the most unstable eigenvalue we found for this 
flow with the same wavenumber pair for the disturbance. The first column of the table gives N , 
the number of radial basis vectors that were used to obtain the eigenvalue given in the next two 
columns, N is also the order of the matrix problem that needs to be solved at each step. 


Most unstable eigenvalue found for the rotating Poiseuille flow (3.12) 
with m = 1, k = —1, Vj = 0.2147, Re = 156.0. 

N 

frequency 

growth rate 

4 

-0.00029 

.00334 

6 

-0.00279 


10 

-0.00284 

.000001 

14 

-0.002847 

.0000001 

18 

-0.002847898 

.0000001379 

22 

-0.002847898 

.0000001378 


The convergence is exponential in N or some power of N and there is no evidence of significant 
roundoff error. The following table lists the corresponding eigenvalue found by solving the 
adjoint viscous stability problem for the same wavenumber pair and baseflow. 


Eigenvalue found by doing the adjoint viscous stability problem for the flow (3.12), 
with m = 1, k = —1, V 1 = 0.2147, Re = 156.0. 

N 

frequency 

growth rate 

4 

-0.00270 

-.000094 

6 

-0.00280 

-.000013 

10 

-0.00284 

-.000004 

14 

-0.002847 

-.0000002 

18 

-0.002847898 

-.0000001379 

22 

-0.002847898 

-.0000001379 


Clearly the agreement between the adjoint and direct results is excellent. Inviscid stability 
results for flows of the form (3.12), obtained using our code also compare well with results in the 
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literature. These results instill confidence in the accuracy of the numerical method and in the 
code that implements it, at least for the case of columnar flows. 

4. Stability Results for the Vortex Breakdown Model Flows. 

In this section we will present the results obtained to date for the stability of the midbubble 
columnar, (1.9,10) and the wavy vortex (1.3,4) flows. Although these flows are inviscid we will 
consider their stability to both viscous and inviscid disturbances (i.e. we will solve the linearized 
Euler and the linearized Navier-Stokes equations for these flows). The justification for doing a 
viscous analysis is that the "real" flow is of course, viscous. Moreover, the inclusion of the higher 
order dissipative terms eliminates certain technical difficulties that arise due to the existence of 
critical layers in the neutrally stable eigenfunctions for columnar flows (Drazin & Reid [1981]). 

We will begin by presenting the viscous results for the midbubble columnar flows. We found 
that thirty Tadial vector modes (N = 30) were adequate to resolve the most unstable eigenmode 
(i.e. the mode whose eigenvalue had the largest imaginary part) to three decimal places for these 
flows at low Reynolds numbers and that this number increased as the Reynolds number grew. 
Frequent checks were carried out on the accuracy of the computed eigenvalues both by increasing 
the order of the expansion and also by computing the adjoint spectrum for the same set of flow 
parameters. The difference between the most unstable eigenvalue as computed by the direct and 
adjoint versions of the code was always less than 1%. 

Having fixed the number of radial expansion vectors in our system the viscous eigenvalues 
for the midbubble flows depend on four parameters, 


a 


a-(m,k,e,Re). 


(4.1) 
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The base columnar flow (e = 0) was found to be stable to all disturbances. It seems that even for 
very small values of e (values for which there is no recirculation zone in the full two dimensional 
flow) the midbubble flows are unstable. This is documented in the following table which gives 
bracketing values for the critical Reynolds number for various values of e. 


Bracketing values for the critical Reynolds number. 
Various values of e and m = — 1. 

e 

Stable for Re 

Unstable for Re 


Stable for all Re 


SB ! ' 

550 

600 


180 

200 


160 

180 


110 

120 


60 

80 

aB 

42 

45 


For large enough values of e disturbances having both positive and negative azimuthal 
wavenumbers can destabilize the midbubble flows with the negative azimuthal modes giving rise 
in general to the largest values for the growth rates. In particular disturbances with azimuthal 
wavernumber, m = — 1 were found to be the most dangerous. This is shown in the following 


table. 
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Bracketing values for the critical Reynolds number. 
Various values of m with e = .03. 

m 

Stable for Re 

Unstable for Re 

-i 

42 

45 

-2 

90 

100 

-3 

700 

800 

-4 

1200 

1400 

-5 

1400 

1600 


For fixed values of m and e, a two parameter ( k , Re) study was carried out. With e = .03 
and m = — 1 we obtain the stability diagram shown in Figure 5. The stability boundary appears 
to be a parabolic curve which is markedly asymmetric with respect to the k = 0 line. Within this 
curve the base flow is unstable to a range of axial wavenumbers; however, there is a ’'tongue’' of 
stable wavenumbers that gradually thins out as the Reynolds number is increased. The k = 1 
mode is the fined one to be excited, this does not happen until Re = 4500 (approximately). 


We now consider the stability of the midbubble flows to inviscid disturbances. The invis- 
cisid stability of columnar flows is governed by an ordinary differential equation, the Howard- 
Gupta [1962] equation. A number of analytic results obtained from this equation exist in the 
literature. Leibovich & Stewartson [1982] showed that a sufficent condition for the instability of a 
columnar flow is that the function, 


F(r) 3 V(r)A'(r)[ A'(r)r(r) + FT( r) J 


(4.2) 


be negative somewhere in the domain of interest. (A is the angular velocity, — V and T is the cir- 


culation rV.) 
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Th function F(r) is easily evaluated for the midbubble profiles and this has been done for a 
range of values of e. It was found that F first becomes negative only when a critical value of e is 
reached, e = 0.132. Consequently the midbubble flows are guaranteed to be unstable by the 
Leibovich-Stewartson criterion for values of the parameter e slightly below those needed to pro- 
duce a recirculation region in the full two dimensional flow (recall this happens for e = 0.155). 

A normal mode stability analysis was carried using the divergence free expansion method 
and the numerically obtained results confirm and somewhat extend the predictions of the 
Leibovich-Stewartson criterion. Once again we found that the disturbances giving rise to the 
largest growth rates had azimuthal wavenumbers, | m| =1. Figure 6 shows how the maximum 
growth rate varies (almost linearly) with the amplitude parameter, e. In the figure we can see 
that the normal mode analysis extends the previously obtained results in that midbubble flows 
with e < 0.0132 are also found to be unstable, even though F(r) > 0 for these flows. The max- 
imum growth rates also increase with the size of the axial wavenumber, | k \ , as shown in Figure 
7. 


We conclude that the midbubble flows are definitely unstable on both viscous and inviscid 
grounds for values of the parameter, e below those needed to produce a reversed flow region in the 
full two dimensional wavy flow. Moreover the most destabilizing disturbances have azimuthal 
wavenumbers, | m| = 1 and short axial period, (| k | > 1). The unstable inviscid eigenfunctions 
tended to have regions of steep gradient near the origin and this made their resolution difficult. 

While these midbubble stability results we have just reported tend to support the conjec- 
tures made about vortex breakdown, they are not encouraging for the numericist seeking to inves- 
tigate the stability of the cnoidal wave flows (1.3,4). As we noted earlier, the stability equation 



- 22 - 


for these flows do not separate in z and consequently we have to solve 0(7^(2/f-(-l)) matrix 
eigenvalue problems for each flow. It would seem to be necessary to include modes having a short 
axial period in our expansion, u^xm which means that K will be large. The inclusion of these 
modes is also dictated by physical considerations. We should like the disturbance to include 
modes that scale with the dimensions of the bubble and in fact visualization studies indicate that 
the asymmetric unstable modes do have short axial periods. However, our experience with the 
midbubble flows show that these unstable modes are difficult to resolve radially, consequently the 
number of radial vector functions, N in our expansion must also be large. 

With these constraints the normal mode analysis of the cnoidal wave flows becomes prohibi- 
tively expensive, (reacall that the number of operations needed to extract all the eigenvalues of a 
matrix is proportional to the cube of its order). To alleviate the cost problems most of the runs 
for these flows were done with the axial period of the cnoidal wave fixed at 1 [H = 0.5, units are 
tube radii). While such short period flows violate the assumptions needed to produce the solu- 
tions (1.3,4) it was hoped that these flows would be unstable to disturbances with smaller values 
for K. Indeed convergence studies indicate that adequate resolution in the axial . direction was 
obtained with K « 10. Up to 40 radial modes were used in the study, leading to an 0(840) 
matrix eigenvalue problem when K = 10. (The calculations were performed on an Floating Point 
Systems 164 series vector processor at Cornell University, using a vectorized version of the QR 
algorithm, optimized for this machine and using some 16 megabytes of memory.) The adjoint 
problem was also solved on each run and only those eigenvalues which agreed well between the 
adjoint and direct calculations were considered. 

Unfortunately the short period flows behave less like the midbubble flows and more like the 
underlying, stable columnar flows. This is indicated in Figure 8 which plots the least stable 
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growth rates found for the various values of e as the Reynolds number is increased (m = —1 in 
this plot). All these short period cnoidal wave flows are stable, though marginally so. The dis- 
turbances with | m| = 1 are again the most unstable. Figure 9 shows the least stable growth rates 
found for some other azimuthal wavenumbers. The inviscid stability runs that were performed 
also failed to turn up any definite evidence of instability for these flows. 

We have considered the viscous and inviscid stability of some axisymmetric flows which are 
said to model the bubble type of vortex breakdown. The investigation was carried out by 
expanding the perturbation velocity in terms of the new set of divergence free vectors presented 
in section 2. The stability results for the midbubble flows support the conjectured mechanism for 
breakdown and it was found that the most dangerous disturbances have a short axial period and 
azimuthal wavenumbers, | m| =1. The two dimensional flows we considered all had rather short 
axial periods and these do not model the physical phenomena well. No conclusions as to the sta- 
bility of these flows can be drawn because the normal mode analysis used here can only prove ins- 
tability (by finding a growing disturbance). No evidence of instability was found but this may 
well be because we failed to include enough axial modes in our expansion for the disturbance. 
The study clearly points out that linear stability investigations for complex base flows are far 
from trivial from a computational point of view. 
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Figure 2. Streamlines (1.3), f = .02, H - z. 
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Axial velocity profiles for the midbubble flow 
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Stability diagram for the midbubble flow with e = .03 and m = — 1. 
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Figure 6. Maximum growth rates plotted against e, m = — 1 and k = 6. 
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Figure 7. Maximum growth rates plotted against axial wavenumber, (inviscid analysis for 
the midbubble flow), m = — 1. 
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